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NUMERICAL METHODS FOR INVISCID/VISCOUS FLUID FLOWS 
T.M. El-Mi stikawy and M,J. Werle 

ABSTRACT 

A study of numerical schemes for solving viscous fluid 
flow problems with sizable regions of predominantly inviscid 
flow is presented. Difficulties associated with the familiar 
central difference approach for such problems are analyzed 
and alternative finite difference approaches employing wind- 
ward concepts are presented. In addition, difference relations 
based on exponential operators are developed. All such 
schemes are demonstrated and evaluated through application to 
the case of Falkner Skan flow with blowing - a problem in 
which a sizable region of predominantly inviscid flow 
developes near the injection surface that traditionally causes 
numerical difficulty. From these studies it is concluded 
that the exponential-box-scheme (EBS) developed here provides 
a definite advantage over the other schemes studied. 
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I. INTRODUCTION 

The s\abject of this study are those numerical methods used 
for solving viscous fluid flow problems with particular interest 
centered on problems where there exists distinct inviscid and 
viscous regions within the calculation domain. For example, 
as depicted in Figure 1, within the supersonic flow field over 
a reentry body where surface ablation occurs , an inviscid region 
is^ formed near the body surface xmder a "blown off" boundary 
layer, which itself is driven by another inviscid region close 
to the bow shock wave. To date, finite difference schemes for 
numerical solution of the governing partial differential equa- 
tions have encountered difficulty with this change in flow 
character apparently due to the diffusion character of their 
difference operators. This difficulty becomes especially 
troublesome when implicit finite difference schemes are employed. 
Although such schemes suffer no integration step size liinitation 
for stability sake, they -are found (see Roache, Ref. 1)- to be 

• ‘ ^ f I -• 

Vt I ' i' 

r’estfrbted by. "accuracy" requirements to virtually the same 
step size ranges as explicit finite difference schemes. Thus 
a need exists to remove this limitation. 

Numerical schemes based on central difference ^concepts 
are found to produce large truncation errors of an oscillatory 
type that do not belong to the physical problem, but rather to 
the difference equations themselves (see Hirsh and Rudy, Ref. 2 
and Roache, Ref. 1) . For want of a better term, Roache (Ref. 1) 
called these oscillations "wiggles" and he gave an explanation 
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for their occurrence based on asymptotic concepts and he 

described a mechanism in which these oscillations occur in a 

simple convection-diffusion problem. The basic difficulty 

appears to be centered in the difference representation of the 

convective process and no completely satisfactory remedy is 

proposed. One solution to this difficulty is based on. the use 

of an active step size control procedure as discussed by Roache 

and used by Liu and Chiu (Ref, 3) . Alternatively, artificial 

viscosity like terms can be injected in the difference 

equation to eliminate these oscillations [e.g. Briley and 

MacDonald, Ref. 4] or one can employ upwind differencing for 

representing the convection terms in the governing equations 

[for examples see the methods used by Garrett, Smith and 

Perkins (Ref. 5) or by Raithby (Ref. 6)]. In addition, it 

has recently been found that numerical schemes based on series 

expansions of the coefficients appearing in the governing 

differential equations result in exponential operators that 

are inherently devoid of this oscillatory difficulty. Such 

* 

schemes were recently proposed by Spalding (Ref. 7) and were 
used by Roscoe (Refs. 10 and 11) and Chien (Ref. 12) for solving 
the two and three dimensional Navier-Stokes equations for a 
simple geometry . 


* A more general approach leading to the use of exponential 
operators was presented by Pruess (Ref. 8) and later by 
Rose (Ref. 9) . 
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In the present work, a detailed study of the problem of 
numerical oscillations and the. different methods used to over- 
come this difficulty is carried out. Finite difference schemes 
based on the addition of artificial viscosity and the use of 
second order accurate windward differences .are presented. In 
addition, schemes based on the method of expanded coefficients 
(exponential schemes) are studied and a new box version of 
same is presented. Application of five difference methods to 
the solution of the Falkner Skan equations with surface 
injection show a clear indication of the superiority of one 
such method (here termed the exponential box scheme) for such 
problems . 
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II. GQVEPimTG EQUATIONS 

1 . General Concepts 

In the numerical solution of the yiscous/inviscid fluid flow 
problems by implicit techniques, the governing equations, which 
in general are partial differential equations , are first reduced 
to two point boundary value problems in a direction perpendicular 
to the principal flow direction. This is, for example, how 
the non self similar two dimensional boundary layer equations 
are solved. Also, an ADI solution of Navier-Stokes equations 
adopts the same treatment. Attention then need initially only 
be focused on the development of numerical techniques for 
ordinary differential equations with split boundary conditions . 

2 . The Falkner-Skan Equations 

A. model for the ordinary differential equations resulting 
from this treatment of the viscid/ in viscid fluid flow problems 
are the mass injection version of Palkner-Skan equations for 
the self similar boundary layer given here as 

+ F = 0 (la) 

- VF^ - B (lb) 

where F is the longitudinal velocity, V is the nomal velocity 
and e is the pressure gradient coefficient. The boundary 
conditions for present interest are given as 

F(0) = 0 and V(0) = V’ (2a) 

w 

with F(n) ^1 as n « ' (2b) 
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The first term in the momentum equation (lb) represents the 
diffusion effect while the second terra represents normal con- 
vection effects. The other two terms represent the longitudinal 
convection and pressure gradient effects. However in this one 
dimensional formulation of the problem they appear as source 
terms . 

I 

In regions where the flow is basically inviscid (such as near 
a wall with large blowing) the diffusion effects are expected to 
disappear and a model of the resulting ordinary differential 
equation representing the momentum balance will be 

- VF^ - = - g (3) 

A numerical scheme that is used for solving fluid flow 
problems where there are transitions from inviscid to viscous 
regions should be able to accommodate this change in the 
character- of the governing equations. 

3. A Model Diffusion-Convection Problem 

A general differential equation that involves both the 
diffusion and convection effects can be put in the form 

F^^ - (a+b)F.^ + abF = c (4a) 

In light of the preceding section, equation (4a) can be looked 
at as a momentum conservation equation where the first term repre- 
sents diffusion effects while the second, noimnal convection 
effects. For initial interest, only normal convection effects 
will be considered by taking b = c = 0 and a = C, a constant. Thus 



C F = 0 
n 


(4b) 
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4, Nximerical Concepts 


There are two methods that can- be --used' ;tc?= .represent ' the dif -7 
ferential equations by difference equations. In the first, the finite 
method. (Ref. 13), difference operators are obtained 
by a Taylor's series expansion approach and these are used to 
replace the differential operators that appear in the governing 
equations . This method is straightforward and provides a formal 
order of accuracy for each such scheme in the limit of small 
step size. In the second, the expanded coefficient method, 
the differential equations are formally integrated over a small 
interval in which its coefficients are approximated by Taylor's 
series representations resulting in a new system of difference 
equations generally involving exponential operators . 


Ill . FINITE DIFFERENCE OPERATORS FOR CONVECTION/DIFFUSION FLOWS 


1. The Central Difference Scheme 

Consider a portion of the flow field depicted in Figure 2 
with j-1, j, and j+1 points of a grid of uniform spacing An, 
and let "m" be the midpoint between the j'-l and jth point and 
"n" 'the midpoint between the j - and g-l^rth point. 


To obtain the central difference representation of equation 

(4b) one can use the Taylor's series expansion to evaluate the 

derivatives F and F at the jth point in the central difference 

n n n 

form ' 



F 


3 +^ 

2An 




F 

nnn j 


+ o(An'^) 


(5) 
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and F 


( 6 ) 






-1 2 
- ^ An F 
12 nnnn. 


+ 0 (An *) 


so that upon evaluating, equation C4b) at the jth point and 

2 

neglecting terms 0(An ) / one obtains the central difference 
equation 

(1 - i An C) - 2F. + (1 + I An C) F,_i = 0 

The exact solution of a difference equation of the form 


«iF-,T+aF +«F = 

1 ]+l ^2 i “3 1-1 


can be written as 


P = + Ej X3 + P 


where and X 2 are solutions of the- characteristic equation 


cc-| X + ^3 ~ ^ 


and R 2 are arbitrary constants determined by the boundary 
conditions and is a particular solution of equation (8) . 
For equation (7) 


= 1 “ ^ An C 
ot 

2 
“3 

“4 
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= - 2 

= 1 + j An c 
= 0 


(7) 

( 8 ) 

(9) 

( 10 ) 


(lla) 

(llb) 

(llc) 
(lid) 



so that 


= 1 (12a) 

1 t .j An 'C 

X = (12b) 

^ 1 - i An C 

and Pj = 0 (12c) 


Thus the exact solution of the central difference equation (7) 
can be written as 


Fj (13) 

where A 2 is given by equation (12b) . 

The exact values of and X 2 can be obtained by comparing 
equation (9) with the exact solution of the differential equation 
(4b) evaluated at the grid points, namely, 

Fj = R^ + R^ (14) 

so that the exact values of the characteristic roots are given as 


X 


1 


1 


X 


2 


e^nC 


Note that in equation (12 b) X 2 will be negative v/hen 


(15a) 

(15b) 


An |cl > 2 ' ( 16 ) 

whereas its exact counterpart, X 2 is always positive. In such 
a case, X^ will change its sign according to whether j is even or 
odd and hence the solution C13) will show an oscillatory behavior 
compared to the monotonic behavior of the exact solution ( 14 ) . 



This oscillatory behavior has been observed in the central' 
difference solution of the Falkner-Skan equations near the outer 
edge of the boundary layer and in many other problems when a 
condition similar to condition (16) and obtained in the same 
way is locally satisfied. 

Thus it is seen that these oscillations belong to the 
difference equations and are a part of their exact solution. 

This eliminates the round off errors due to the lack of 
diagonal dominance from being a cause of these oscillations. 

Roache ('Ref. 1) attributes these oscillations (wiggles) to a singular 
behavior of the difference equation (7) when the cell Reynolds 
number An[cl becomes large (An|c[ >2) compared to the 
singular behavior of the differential equation (4b) in the limit 
when the effective Reynolds number |cj becomes infinitely large. 

To give a physical interpretation to this "unstable" 
behavior of the difference operator, it is convenient to employ 
a control volume (or conservation) representation of these 
equations. The central difference equations (7) can be obtained 
through a control volume approach by first integrating equation 
(4b) over the control volume [m,n] to get 



(17), 

(18a) 

(18b) 


<3 



(18c) 




+ o(An^) 


F 

n 



+ O(An^) 


(18d) 


Note that in this representation of the convection terms there 
is an explicit downstream contribution to the momentum/ carried 
by the flow as it moves with the velocity C. This downstream 
contribution is a characteristic of a central difference repre- 
sentation of the convection effects that is expected to cause 
large errors and instabilities of the scheme. 

In fact, the oscillatory behavior described in the previous 
section is how these instabilities manifest themselves. Note 
that according to condition (16), the oscillatory behavior 
is encountered when the step size An is large in which case the 
downstream contribution to the convection effects comes from far 
points, or when the normal velocity C is large in which case 
the downstream contribution itself becomes large. 

Condition (16) for the oscillatory behavior suggests the 
control of the step size so that 

An |c| £ 2 (19) 


for no oscillations. 

For the Falkner Skan equations where the normal velocity V 
varies from point to point a condition equivalent to condition 
(19) should be locally satisfied. Thus, in this case the control 
of the step size can take one of the following forms : 
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1. - using, with a uniform grid, a step size that uniformly 
satisfies the condition for no oscillations. ' Thus even in regions 
where larger step sizes can be used, the step size should be less 
than or equal to that corresponding to the severest case. This 
means a larger n\imber of mesh points and consequently greater 
computational time and storage requirements . 

2. using a variable step size that locally satisfies the 
condition for no oscillations as used by Liu' and Chiu (Ref. 3) for 
massive blowing problems. Actually use of Keller's box scheme 
(see Blottner, Ref. 14) , which is a version of the central 
difference approach is directly adaptive to this method. How- 
ever such an approach is unsuitable to nonlinear equations since 
in an iterative procedure the solution changes from one iteration 
to another. 

Two other approaches can be used. The first is the use of 
non-centered differencing for the convection terms which will 
be discussed in a later section. The second is the addition of 
explicit artificial viscosity terms to the governing equations . 


2 . Artificial Viscosity - Corrected Central Difference Scheme 

In the second method terras of the form vF are added to 

nn 

the differential equation, e.g. equation (4b) becomes 


(1+vO F - CF = 0 

rm Ti 

with the central difference representation given by 

F 


( 20 ) 


( 1+v ) 


- 2F. 
1+1 ~l 


'±± - 


C. 


An' 


P — P 

1;^ :.i - i, 

2An 


= 0 


( 21 ) 
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whose characteristic roots are 


= 1 - (22a) 

1 + V + i An C 
^ ^ 

^ l + V'-i-AnC (22b) 

J 

So, the condition for no oscillations will be 


An [cl £ 2(l+v) 


(23) 


The artificial viscosity v is usually constructed to be of higher 
order so that the accuracy of the method remains second order, 

A more consistant method to introduce this artificial vis- 
cosity, can be based on a second order correction to the central 
difference representation of the derivatives. In equations (5) 

and (6) one can relate F and F to F through 

nnnj nnnHj nrij 

differentiation of equation (4b) to get upon substitution in 
equation (4b) evaluated at the jth point 


p 2 f + p 

(1 + ^ '^'>2 c^) 

An^ 


- c 


P - F . 

liri 

2An. 


0 (24) 


Serice , comparing with equation C21). the correction introduced 
is equivalent to an artificial viscosity 


1 2 2 
V = A.2 c" 


which is second order in An and hence does not harm the order 
of accuracy of the scheme. 

The difference equation (24) can be put in the form 


[1 - i C + ^ An^ ^^^'^j+1 ^ An^ C^JPj 

+ [1 + I An c t ~ An^ C^]F. ^ = 0 


(25) 
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which has the characteristic roots 


1 2 12 2 
(1 + I An C) + 

(1 - I An C)2 + ^ An^ 


(26a) 


(26b) 


Thus X 2 always positive and the solution is non oscillatory. 

From equations (24) and (26b) , it is seen that the smallest 

2 2 

amount of artificial viscosity proportional to An C that can be 

1 2 2 

added to prevent the oscillations is An C which can be 

looked at as that part of the truncation error that is responsible 

for creating the oscillations. This suggests that the L.H.S. 

of equation (24) can be used to represent the differential 

expression F - CF whenever it appears in a differential 
Tin j n j 

equation instead of introducing a complete correction to the 
difference equation which requires differentiation. 


3 . Upwind Schemes 

The analysis of the previous sections has shown how a 
central difference representation of the convection terms can 
produce the problem of the oscillatory behavior of the numerical - 
solution, and has related that to the explicit downstream 
contribution to the convection effects . 

The idea behind upwinding is to correctly represent the 
convection terms so that the information carried by the flow as 
it moves is the upstream information only. 
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A first order accurate upwind representation of equation 
(4b) can be obtained if in equation (17) the diffusion terms are 
represented by equations (18a/ b) and the convection terms are 
represented by 


and 


F = F. T + 0(An) 
m 


F = F . + 0(Ati) 
n 3 


(27a) 

(27b) 


when C > 0, i.e. the flow is upward, and by 


and 


F = F . + 0(An) 
m J 

F = F. ,, + 0(An) 
n 


when C < 0, i.e, the flow is downward, 
equation can be put in the form 


(28a) 
(2 8b) 

The resulting difference 


[1- I An (c-jcj ) -[2+An|c|]Fj +[1+ j An(C+ jcDjFj^^ = 0 


(29) 


which can be also obtained through the Taylor's series expansion 

approach if equation (4b) is evaluated at the jth point and a 

central difference representation of F and the following 

nn^ 

upwind difference representations of F are used 


F. - P. , 
3 1-1 

An 


+ 0(An) 


n. 


if c > 0 


(30a) 


F . - F . 

+ o(An) if C < 0 


(30b) 
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The 

characteristic 

roots of equation (29) are 





(31a) 

and 

1 

\ — 

+ i An (C+ |ct ) 

(31b) 

^2 

^ 1 

- i An(C-lc ) 


In both cases C > 0 and C < 0, ^2 positive and the solution 

is nonoscillatory which supports-, the upwinding idea. 

Applied to the Palkner Skan momentum equation, where the 

normal velocity V changes from point to point and may even change 

its direction, the first order accurate upwind scheme described 

above will have to follow the direction of V and to switch 

when its direction is changed so that the information carried by 

the flow is always the upstream information. 

Another approach to obtain a first order accurate upwind 

representation of equation (4b) is to evaluate this, equation at 

point "m" when C > 0 and to use the following second order 

accurate representation of F and first order accurate repre- 

m 

sentation of F 

m 

P . - F . 

P = -J— hk + o{An^) (32a) 


F 

nn 


m 





+ o(An) 


(32b) 


when C < 0 equation (.4b) is evaluated at point "n" and the 
following representation of F and F are used 
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F 


n 


n 


F 

nn 

n 


j+1 


- P-. 


-An--- 


+ o (-An ) 


^ 1+1 ^ ^ ■ 


+. 0-(An ) 


•^(33a). ■ 


(33b-)- 


in the case of equation (,4b) this- approach, will lead to. the same ' 
difference equation (29) asthe previous approach. However. this 
■will not be true if C is variable dr in the "case of a nonlinear 
equation ‘ like the .Falkner. Skan. momentum equation. In fa'ct, 

•there are two main differences between the -previous .twp approaches. 
The first difference ,:xs in, the order of accuracy in which, the 
diffusion and' convection -terms are' represented. For the first 
approach, the diff-usion terms are second order accurate while 
the convection terms are first order accurate. For the second 
approach --the diffusion ‘terms are first order -accurate whi.le 
the. convection .terms are second order accurate. . Thus .in a region 
of inviscid flow where", the. ..diffusion term disappea.rs from- the ‘ 
momen-tum equation the"' second approach will lead- to a second 
order accurate scheme- while’ 'the first approach will lead to a- 
firs-t order accurate scheme.' The second difference is in the 
control vol-ume for which the difference equation is the momentum 
balance equation. For the first approach the control volume is 
the interval between the -two midpoints "m” and "h". For -the 
second approach it 'is the interval between j-1 and j if C > 0 
and the interval between j, and j+1 if C < 0.. 

A special treatment at the switch region should take place 
to preyent excessive errors due to the nonconservative nature 
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of these two upwind schemes. This is especially required for the 
second approach where the momentum conservation for a control 
volume near the switch point is not considered. 

Corrections to the first order upwind schemes described 
above to be second order accurate can be introduced by considering 
the next term in the Taylor's series expansion representation for 
the first order accurate derivative involved. These corrections 


can be introduced either in an explicit sense using values from 

the previous iteration or in an implicit sense . 

In the discussion that follows, only the case when C > 0 

will be considered. Similar treatment can be applied to the case 

when C < 0 . Consider first the windward approach that corrects 

the first order representation of the windward convection 

operators given by equations (30a) and (30b) . For this approach 

note that the first derivative F can be expressed as 

■ ■ — 

F = F +4atiF + 0(An ) (34) 

3 m m 

In equation (30a) , only the first term was considered' and was 
‘represented with second order accuracy. To get a second order 


representation of F^ , the second term should also be considered. 

However, since in that term F is multiplied by An one can, 

with first order error, evaluate F at any close point and 

nn ^ 

still have second order accurate representation of F 

nnj 

One choice is to evaluate F^^ at the (j-l)th point. This leads to 



(35) 
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where central differences have been used to represent 

{ j-1) 

and the "barred" values refer to the previous iteration or 
station values (see Atias , Wolfshtein and Israeli, Ref. 15). 

This method is consistent with the upwind idea in that 
upstream information is used to evaluate the correction, however, 
a special treatment near the end points will be needed. 

A second method is to evaluate F at the jth point v/hich 

Tin 

leads to 


F 

ri 


j 


- F 

TrT 




1 . 
2 


F 


j tl 


- 2P. + F 


2ZI 


o(An^) 


(36) 


An' 


Eq. (36) was obtained by Khosla and Rubin (Ref. 16) by rearranging the 

central difference representation of F^ [equation (5) ] . A 

scheme based on this method will have the same magnitude and type 

of error as the central difference scheme within the converged 

solution since they are basically identical in representing the 

derivatives F and F . Thus this method will give an oscillatory 

n . nn • 

3 3 

solution if the central -difference solution is oscillatory. 

However this scheme has the advantage over the central difference 
scheme that it is always diagonally dominant (see Appendix A) . 

A third method is to use the governing equation to represent 
F in terms of the other terms appearing in the differential 

equation. From equation (4b) one can write 


nn. 


= C F 


n 


m m 

Substituting into equation (34) one gets 


(1 + i An c" ) F + O(An^) 


18 



In difference form 


n. 


(1 + ^ An C) 


F. - F. T 

_J irl 

An 


+ 0{An ) 


(37) 


which leads to the difference equation 
1 


- 2F. + P. 
3+1 3 3-1 


1 + J An C 


An' 


F. - P. , 

_3 Hi 

An 


= 0 


(38) 


This is termed the grid point upwind scheme. 

Note that, with this correction, the so-called implicit artificial 

viscosity associated with the first order upwind scheme has been 

removed by reducing the effective viscosity in the difference 

equation. This suggests that the L.H.S. of equation (38) can be 

used to represent the differential expression F - CF when 

Tin j n^ 

C > 0 whenever it appears in a differential equation. This 
can be shown to be second order accurate. 

Considering the two cases C > 0 and C < 0 a combined 
difference equation for equation (4b) can be put in the form: 

[1 - An(C-lc[)+ i An^(C-lcl)^]Fj^^ -[2+Anl,c! + ~ An^C^lF^ 

+ Cl + I An(C+]cj)+ I An^(C+|cl)^]Fj_^ = 0 (39) 

whose characteristic roots are 


- 1 (40a) 

1 + y An(C+[ci)+ y An^(C+|c|)^ 

A„ = 1 f (40b) 

1 - y An (C-|cl)+ ^ An (c-jcj)^ 

^2 it equation (40b) is always positive. Moreover if C > 0 
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then X 2 > 1. While if C < 0 then ^2 ^ 1 which is in accordance 
with the behavior of in equation (15b) . 

Consider now the second windward approach that evaluated 
the goveiTiing equation at point “m" (or ”n") using a second 
order accurate representation of the convection term and a 
first order accurate representation of the diffusion term. 

To correct this approach to second order note that the second 
derivative F is given, to second order accuracy by 




~ An + O(An^) 

2 Tinn^ 


(41) 


where, again, a first order accurate representation of 
will still lead to a second order accurate representation of 



Taking 


F 

nnn j 


nnn 


m 


0 (An) 


and using the derivative of the differential equation evaluated 
at point "m" to evaluate F one gets , in the case of equation 
(4b) , 


nri. 


m 


1 + i AnC 


F + 0(An ) 


With central difference representation of F this becomes 


nn. 


nn 


m 


1 + J AnC 


F. , - 2F. + F._. , 

^ + o(An^) 
An 


(42) 
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and the difference equation (.termed the midpoint upwind scheme) 
will be 




- 2F. 

1 


^-1 _ 


F. - 



1 + 


AnC 


An' 


An 


F. T 


(43) 


which is the same as equation (38) and hence equations (39) and 
(40) will apply. However, this is not true for more general 
differential equations than (4b); i.e. the difference equations 
will in general, be different. 

Note that a family of three point finite difference equations 
can be obtained by applying the same procedure described above 
but for a general point rather than the midpoints (see Appendix B 
for details) . 


4 . Summary and Comparison of Finite Difference Schemes 

An assessment of the mmerical schemes described above can 
be made by comparing the characteristic root ^2 for each scheme 
with the exact value X 2 as given by equation (15b) . 

In Figure (3) , for the central difference (C) , corrected 
central difference (CC), and the second order upwind (UW) schemes 
as given by equations (12b) , (2Sb) , and (40b) respectively and 
X 2 of the exact solution are presented. 

It is clear that for small An|cj all schemes have reasonable 
agreement with the exact solution. This agreement is best in 
the case of the corrected central difference scheme which is 
fourth order accurate for this special problem. 

For Arilcj >2 the central difference scheme has negative 
^2 which produces the oscillatory solution and it has a singularity 
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at AnC =2. ^2 for the corrected central difference scheme 

is always positive however when AnC = i/12* reaches an 

extremum after which it asymptotically approaches ±1 to ±“. 

Thus this scheme is expected to produce large errors for 
Anlcj > /l2 . 

The two branches of the upwind scheme (for positive and 
negative C) each follows the behavior of the exact X 2 in its 
region of application. On the other hand, out of their 
region of application [designated by the downwind scheme (DW) ] 
jthey encounter large errors, and would not be reliable. 

IV. THE METHOD OF EXPANDED COEFFICIENTS FOR 
CONVSCTION/DIFFUSION FLOWS 

1. The General Concept 

Consider the general second order differential equation 
given by equation (4a) as 

F - (a+b) F + abF = c (44) 

HT] n 

where a, b, and c are in general functions of n and P. Instead 
of using Taylor's series expansions to obtain finite difference 
representations for F, F , and F , one can instead expand the 

T1 TIT) 

coefficients a, b, and c about a point in the interval under 
consideration. The first terms in these expansions are the 
values of those coefficients at that point. The differential 
equation that results if only these 'first terms are kept will 
be an ordinary differential equation with constant coefficient 
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that can be solved analytically. One can use this exact solution 
to obtain finite difference equations as will be described in 
the next sections. 


2. The Three-Point Exponential Scheme 

If a, b, and c are expanded about point j and the resulting 
differential equation is assumed to apply in the interval 
[j“l/ j+1 ] t an exact solution of the resulting constant coeffi- 
cient differential equation will be given as 


a .Ti 


b .n 


F = Ae^ +Be^ 


a .b . 
3 3 


(45) 


Using this solution to evaluate F . , , F . , and F . , , one obtains 

^ j“l j j+1 

three equations that involve the two arbitrary constants A 

and B which can be eliminated to give the following difference 

equation : 


a . An b . An a . An b . An 

F.^, -(e^ )F.+e^ e^ F., 

3+1 3 3-1 


a .b . 
3 3 


{1 -(e 


a^An b.An a. An b.An 

+ e^ )+e^ e^ } 


(46) 
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A Taylor's series expansion of and S'j_j_ about point j and 

the expansions of the exponentials involved for small An shows 
that this difference equation is a second order accurate 
representation of equation (44) ,, (see Appendix C) . 

For the special case when b = c = 0 and a = C , 
equation (44) reduces to equation (4b) and the difference equation 
(46) reduces to 

CAn 

Fj +1 - (1 + e ) Fj + e = 0 (47) 

which has the exact characteristic roots 




and ^2 = e 


CAn 


Roscoe (Ref. -100 used equation C47) to derive 
representation" in the form 


C Tt:, /I , 

-1) ® J 


a "unified difference 


(48) 


to represent the differential expression F -CF with second 

nn n 

order accuracy (here teimied the RUDR) . 

If the exponentials in expression (48) are expanded and first 
order terms in An are neglected this expression reduces to the 
first order upwind expression for positive C. The second order 
upwind expression follows if instead second order terms in An 
are neglected. To get the corresponding expressions for negative 

f f t 

C one first multiplies both the numerator and the denominator 
of (48) by e and then expands' the exponentials. 
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The central difference expression follows if the exponentials 

~Ari C /2 

resulting after multiplying and dividing by e are 

expanded to second order. 

Instead of using (48) to represent one 

can use equation (46) which gives a consistent representation of 
the whole differential equation at the expense of the complications 
of evaluating the auxiliary roots a and b. 

Consider now the limiting case when a with b and c/a 

fixed. Equation (44) reduces to 

-a F + abF = c (49) 

n 

which corresponds to an inviscid flow of velocity a. In this 
case, equation (46) becomes 

b . An c . b . An 

-F. + e ^ P. T = — j— (e ^ -1) (50) 

1 1-1 a .b . 

j j ^ ^ 

which also follows if the approach used to derive equation (46) from 

equation (44) is applied to equation (49) in the interval [j-l)j3 

with the coefficients a, b, and c evaluated at point j. This 

scheme is first order accurate (see Appendix C) . An equation 

that involves F. and F. , similar to eguation (50) follows if 

3 3+1 

the limit as a -»■ -« is considered instead. Thus if this 
scheme is applied to a problem where there is a transition 
from viscous region to an inviscid one its accuracy will change 
from second order to first. 
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3. The Exponential Box Scheme ~ 

To overcome this difficulty let the coefficients in equation 
(44) be evaluated at the midpoint ”m" and "the resulting differ- 
ential equation be valid in the interval [j-1, j]. Then the 
exact solution can be obtained in the fom 

c. 


a.| n b-, ri 

F = e e H- ^ 


(51) 


1*^1 


where the. s\±iscript "1" -refers to the interval tj-l/jl* 

With ri=0 and n==An at the (j-l)th and jth points respectively 
one gets the following expressions for F , F. and P 

J 


F . T = 
3-1 

^1 

«U 

®1 

4- 

‘^l 

(52) 





a. An 


b. An 
e 

4. 

=1 

(53) 

r . 

3 

^1 ® 


1 


P 

Three other 

aiAn 

equations 

b. An 

+ b^B^e 

for. the interval [j 

(54) 

j,j+l] can be obtained 

in the foinn 







F . 
J 

^2 


®2 • 

4. 

^^2 

(55) 


r 

^2^2 

p = 

^ j+1 

a,An 

a ^ 


b.~An 

B2 e 

4* 

°2 

(56) 




^2^2 

P 

■n 

a-,A., 

+ 

b,B„ 



(57) 
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where now the subscript "2" refers to the interval [j, j+1] and 


3 .^, and are the values of a, b, and c at the midpoint "n;" . 


The requirement that F of equations (54) and (57) should be 
the same together with the other four equations make it possible 
to eliminate the four arbitrary constants A^, A 2 and B 2 
to get the difference equation 
-a2An 


(a2*:b2)e 


-a2An 

a-,e e -b« e 

F. -- . - ■ + 


-a^^ATi b^Ari 
e 


-a2An b2An j+i ^2^^^ 

1-e e 1-e e 


} 

-a^Art j 

1-e e 


(ai“bi)e 


bj^An 


(a2"b2)e 


-a2An 


-a,Ari biAri j-l 
1-e e 


-a2An b2An ^ 2^2 
1-e e 


- { 


^2® 


-a2Ari b2An b^Ari 

® -h °2 ^ ® 


1-e 


-a2Ari b2AT) 


^2^2 


-a^An 
1-e e 


aibi' 




bj_An 


-a. An An 
1-e e 


(58) 


If the limiting case of equation (49) is considered/ equation 
(5 8) reduces to 


b^ An b. An 




-F . + e p. - = .--r; ■ (e -1) 

for 

a 

(59a) 


j j-l a^b^ 


b^An bpAn 




and -F. +e F . = _ , (e -1) 

for 

a -J- -CO 

(59b) 

j+1 3 ^ 2^2 





which are second order accurate (see Appendix C) . 
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Thus the scheme represented by equation (58) is second 
order accurate in both the inviscid and viscous regions and it 
switches smoothly from one region to the other following the 
upwind rules as seen from equations (59a and b) . 

In addition to this advantage this scheme has, like the 
previous one, the advantages that it is diagonally dominant for 
negative ab (which corresponds to favorable pressure gradient) 
and that the error terms are independent of the velocity (a+b) 
compared to the other finite difference schemes described 
previously. In fact the error is due to approximating the 
coefficients a, b, and c by constant values and hence it depends 
on how those coefficients behave with n • 
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V. APPLICATIOH. TO EALKNER SKAN FLOWS' WITH INJECTION 


1. General Statement 

Numerical solutions to Falkner-Skan equations with injection 
bomdary conditions given in equations Cl) and (2) were obtained 
using the central difference, corrected central difference, 
midpoint upwind, Roscoe's UDR, the- three- point exponent icil', and 
the- exponential box- schemes • 

The boundary condition (2b) was applied at a sufficiently 
large but finite distance n^* A grid of uniform mesh An and 
N+1 grid points (n^ = W An) was used to divide the region between 
the surface (n = 0) and the boundary layer edge (n = n^) into 
N intervals . 

Equations (1) were solved in an uncoupled sense. The 
momentum equation (lb) was linearized by evaluating the 
coefficients of F^ and F using the previous iteration values 
to give 

F -VF-eFF=-e (60) 

nn n n 

where the "barred" quantities are the previous iteration values. 

Each numerical scheme was used to reduce equation (60) 
to a difference equation of the general form 

“ll ^j-H “2j^j + “3j^j-l = “4j <«1) 

for j = 2, 3, ..., N. The resulting N-1 equations together with 
the two boundary conditions 

= 0 and ^N-fl ~ ^ 
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form a tridiagonal set of N+1 equations that were solved with 
the standard Thomas algorithm. to give F^, j=l, . . N+1. 

Then the continuity equation (la) with the boundary condition 
(2a) was solved for V ^ , j=l, .../ N+1. For the first five of the 
schemes considered (i-o. central/ corrected central/ windward/ 
Roscoe' S' UDR/ TP- exponential schemes') the. continuity equation was 
integrated by the trapezoidal’ rule which results in' the difference 
equation 

An 


= ''j-i - r 


(63) 


For the exponential box scheme the expression for F in each 
interval was introduced into -the continuity equation which was 
integrated analytically to obtain a difference equation for V. 

This procedure was repeatedly applied to relax a guessed 
solution to a converged one. The initial guessed solution 
that was used in all cases was 



(■6 4)' 


Fj = 1 j = 2, . . . / N+1 (65) 

and the corresponding V^'s were obtained by integration of 
equation (la) with the trapezoidal rule as expressed in 
equation (63) . The solution process was considered to be 
converged when the maximum difference between any two corre- 
sponding F's in two consecutive iterations was less than 10~^.* 


* iUi additional condition on V is used. The maximum difference 

between any two corresponding Vs is normalized w.r.t. the old 

value of V if it is greater than unity and w.r.t. unity otherwise, 

and this normalized difference is required to be less than 10"5. 

This condition assures that V has converged to the fifth 'decimal 

place if it is less than 1 and to five significant *digils otherwise. 

■< <1 
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2 . Difference Equations 

(a) The Central Difference Scheme 


Introducing the central difference expressions (5) and 

(6) for F and F into equation (60) evaluated at point j 
rij Tinj 

straightforwardly gives the central difference equation in 
the form 

(1 - I AnVj)Fj_i_^ - (2 +An^ ^ k 

( 66 ) 

(b) The Corrected Central Difference Scheme 

The differential expression F^^ - F^ in equation (60) 

evaluated at point j is replaced by the difference expression 


F, -2F.+F. F -F 

1 * 2 - 2 , ^ j +1 ^ 1-1 _ 1+1 1-1 


(1 t ^ An Vj) 


An' 


V. 

J 


2 A n 


to give the difference equation 


(X - I 4, V. + ij 


-(2+ ~ An^ + An^e F.)F. 
6 J 3 3 


+ a + i An V. + An^ v2)f._ 3_ = -An^s 


(67) 


(c) The Midpoint Upwind Scheme 

In this scheme the difference equations are written at 

N-l midpoints only. The midpoint at which V first becomes 

negative is not included. If at the midpoint "m" (Fig. 2) 

V is positive, equation (60) will be evaluated at "m" to 
m 

give 


F . - V F - gF F 

nri^ ^ '^m ni- m 


g 


( 68 ) 
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Fjjj and are represented by 


run- 


F . + F . 

F = + O(An^) 

in z 


(69a) 


F . - F . 

and F = ^ + O(An^) 

m 


(69b) 


To find the representation of F one uses the relation 


F =F -^AriF +0 (An^) 
nil nri. 2 nnn 

m 3 m- 


(70a) 


with F given by the derivative of equation (60) evaluated 
in 

at "m" as 


nrni 


m 


VP +V F +26FF 
m Tin n n m n 

m m- m m 


Using the continuity equation this becomes 

F = V F - (1 - 2g)F F 

TjTiTi_ m T]n m n 

mm m 

Evaluating F^ and from the previous iteration values and 
substituting into (70a) one gets 

(1 + i An V )F = F + \ An(l-2g)F F + O(An^) (70b) 
z m n n^ n n • 2 m n„ 

m 3 m 


With a central difference representation of F combined with 


nn. 


equation (69a) for F , this becomes 


m 


(1 + - An V )F 
2 m nn 


F . , T - 2F . + F . , , 


m 


An' 


P .-F 

. y An (1-2B)F — ■ 

^ m ^ n 


-f 0{An ) 


(70c) 
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Substituting (69a) , (69b) and (70c) into C 68 ) the following 

difference equation is obtained 


j +1 


- [2 + A’n V, 


m 


1 . 2 -2 
2 ''m 


- i 4n^a-2B)F 

m 


+ An (1+ 


1 

2 


AnV ) BF ]F. 
m m 3 


+ [1+ AnV + i An^ ^ An^(l-2B)F - An^dt i AnV )B F )F. ^ 

m- 2 m- 2 m 2 m m 3 -I 

= -An^B (1 + i An V^) (71) 

2 m 

When V^ is negative, equation (60) is evaluated at "n", the same 
procedure outlined above is applied and an equation similar to 
(71) is obtained. 

In equation (71) , taken to be the average of F^ and 

F. , with second order error. A similar treatment of V was 
D -1 m 

found to cause large errors so that a fourth order accurate 
value of V^ was used. Using Taylor's series expansion, one 
can write 

V. + V. , , , _ . 

V = - 1 -^ — InL - ^ V to(An') 

m 2 0 nn„, 

m 

The differentiation of the continuity equation gives 


V 

tin 


m 


F 

n 


m 


so that 


V 

m 


V. + V. , ^ P. - 

-J izi + 1 An^ 

2 ^ 8 An 


+ 0 (An ') 


( 72 ) 


where the difference representation in (69b) is applied to 

F . 




IS 
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(d) Roscoe's UDR Scheme 


The difference expression (48) with C replaced by Vj is 

used to replace the differential expression F . - V.F . in 

nnD J nj 

equation (60) evaluated at the jth point. The resulting dif- 
ference equation will be 


V. 

1 

v.An 

An (e ^ “D 



-C1+ 



V An 

e )F. 

: 




V.An 

e 


F 


j 



(73) 


(e) Three Point Exponential 'Scheme 

The coefficients in equation (60) are evaluated at the 
jth point so that 


F -V.F -gF.F=- 
nn J n o 

which now corresponds to the 

(44) . The auxiliary roots, a 


3 (74a) 

general form given in equation 
j and bj , are given by 


Sj = Vj/2 + /(Vj/2)^ + gFj 


■ (74b) 


and bj = V^/2 - /(Vj/2)^ + gF^ 


(74c) 


with c j = - 6 . 

The resulting difference equation is then given by 
equation (46). 

Considerations must be taken to avoid dividing by zero in 
the right hand side of equation (46) when one or both of a^ and 
bj are identically zero. To avoid this .difficulty, the R.H.S. 
of equation (46) can be written as 
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R.H.S. = 


a . A n b . A'Ti 

'I' -■ e 1 e'-’^ 


c . 

b . • 3 

3 


which has the obvious limit ■ for b . -j- '0 

3 

b:. An 

, lim r ■- e ■-* 

*J*0 

3 . 


so that 


R.H.S. = - An 


a . An 
■1 e ^ 


If both a^ and b^ are zeros then a limiting process on a^ 
similar to that in (76) puts (77) in the form 

R.H.S, = An^ Cj (78) 

It 'Should be noted that this particular limit treatment 
results in a difference equation- which corresponds exactly to 

J 

the limit form of the differential equation under consideration. 
For example, if b^ = 0, the difference, equation will be 


3+1 


a. An a.An 

(e ^ +l)Fj + e ^ 


a . An 

3 ^ 


which will be obtained if the exponential approach is applied to\ 
the differential equation 

F - a. F = c. (8 0) 

■ nn 3 n 3 

(f ) Exponential Box Scheme 

In the difference equation (58) the roots a^ and b^ are the 
auxiliary roots of the differential equation (60) given in equa- 
tions (74b&c). with the j subscript being replaced by "m" of 
Figure 2 . Likewise the and are the roots corresponding 
to the midpoint "n". For the present application, also 


=1 = =2 =• - « 



Again, special treatment is necessary for cases when any 
of these auxiliary roots vanish. ■ For example,' the coefficient 
of in the right hand side -of equation (58) can be put in the 


form 


a„An b_Ari 
1 - e e 


a An b^An 9 ^ ^ 

1 - e . 1 - i 2 , 

[e + 1 


(81) 


and a limiting process: similar to that in (76) can be applied 

_a,An ^iAti 

if necessary with 1/1 ~ e e dropped from both sides 

of equation (58) if all a^’s and are zeros (i=l,2) . 

As mentioned before, the integration of the continuity 
equation (la) was performed in a way consistent with the 
treatment of the momentum equation. The expression for F in 
an interval [j-1, j] given by (51) is substituted into (la) 
which is then integrated to give 

A^ a-, n B-, b-, n c. 


V = - — e 


'^1'’ 

-^1 


— ri + D-, 


(82) 


where A^^ and are obtained in terms of ^.nd F^ from 

tiii 

equations (52) and (53). Evaluating (82) at the j-1 and 

+* Vi 

j points where n=0 and An respectively and then subtracting, 
a difference equation for V is obtained in the form 

. A- a. An . B, b. An c. 


1.. a^ An . B, b. An 

V. =-V..+ — (1-e )+^(l-e^ ) - 

3 J-1 a^ ^1 


An 


(83) 

To evaluate the midpoint values of F and V which are necessary 
for evaluating the auxiliary roots a^ and b^, the expressions 
in (51) and (82) are used. 
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To avoid the difficulties that arise due to this need of 
evaluating F and V at the midpoints another version of the 
EB scheme is obtained and described in Appendix D. Also, a 
new inversion algorithm that can be used with that scheme 
is presented in that Appendix. 
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VI. NUMERICAL RESULTS AND DISCUSSION 


1. General Statenient 

Three cases of the> Falkner-Skan equations with different 
values of the pressure gradient coefficient 0 and the injection 
velocity were solved over a range of mesh step sizes-- .A 
step size study for the calculated velocity F at chosen points 
was performed to answer two questions: 

1. How long would each scheme considered stay second order 
accmrate as the step size An increases?, and 

2 . How accurate is each scheme compared to the others in 
.the different problems conside-red? 

In addition, attention was focused on the numerical behavior 
of the velocity profiles to determine if local oscillatory 
behavior ("wiggles") was encountered. 

In all cases, the calculated P at the chosen points was 
plotted versus the square of the step size and the resulting 
curves extrapolated to zero step size to provide an "exact" 
value as a basis of comparison. In general, these plots should 
be straight lines for small step sizes since the schemes considered 
were all second order accurate ‘in An- The- answer, to the first 
question above -can therefore be determined by obseifving how long 
these plots remain straight as the step size was increased , 
while the answer to the second question is obtained by comparing 
the slopes of these terminal line segments . 
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2 . 


Case I - The Blasius Squa-fc'ions 


As -a standard of comparison, attention is first directed 

to the familiar flat plate case corresponding to zero pressure 

gradient (B ~ 0) and no injection (V = 0) . The momentum 

w 

equation (lb) reduces to 

F - VF = 0 (84) 

nn n 

The F-velocity profile is shown in Figure 4, where it is seen 
that the thin boundary layer is confined to the region r\ < 6 . 

Figure 5 shows the two individual terms of equation (84). 
demonstrating the complete balance of these terms at all 
points across the viscous layer. 

. The two points n = 1 and, r\ — 2 were chosen to perform 

the. step size study. The first point is a typical point in 

the boundary layer, while the second is the point where the 
diffusion and convection effects have maximum values (see 
Figure 5) . 

In Figures 6 and 7, four step sizes An = 0.2, 0.25, 0.3, 

and 0.5, are used to plot F at n = 1 with straight segments 

joining the points so that any change in the slope can be 
observed. In Figure 6, a straight line is obtained by the 
central difference scheme (C) , while the upwind Scheme (UW) shows 
deviation from the straight line as the step size increases. The 
corrected central difference scheme (CC) which is seen to be more 
accurate than the other two, shows a slight curvature. In Figure 7, 
the three point exponential scheme (TPE) and the Roscoe's UDR scheme 

original PAGB^ 

39 POOR (JUALTIY 


(RUDR)‘ which are identical for this special problem^ give a 
linear curve over the chosen range of An . This curve is very 
close to that of the CC-scheme and in fact represents a straightened 
version of the scheme. This observation (for the RUDR scheme only) 
is found to be true in all the. results obtained in this study. 

The exponential box scheme (EB) is also seen to follow a 
straight line over this range of An and is, slightly less accurate 
than the RUDR and TPE schemes. 

In Figures 8 and .9, one more step size An = 0.4 is used to 
plot F at n = 2, Again, in Figure 8 the CC-scheme shows better 
accuracy than both the G -scheme (which also loses linearity) 
and the UW-scheme (which displays a subtle nonsystematic behavior 
in that its slope increases and decreases from one segment to 
another) . These behaviors of the C. and OW-schemes are presumed • 
to be present for F(l) as shown in Figure 6, although they were 
not observed- within the accuracy of the plots . As shown in 
Figure 9 , all the exponential schemes produce- straight line 
variations over the entire range of An for this point. 

In Tables 1 and 2 , the detailed velocity profiles are 
presented for a step size of .0.5, where local oscillations are 
observed in the central difference calculations near the outer 
edge when jv| > 4, i.e. when An [v| > 2 which is in accordance 
with condition (16) . 
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3. -Case II - Moderate Injection at an Axisyitimetri'c Stagnation Point 


Here an injection velocity of = 4 is used with a pressure 
gradient coefficient of 3 — 0.5/ so that a region of primarily 
inviscid flow near the surface is created. In this inviscid 
region (n < 5) , the velocity profile is shown in Figure 10 to be 
nearly linear. Above this region (n > 5) , the diffusion term 
begins to play an equal role in the force balance represented 
by the momentum equation. This classification is more clear 
in Figure 11 which presents the. individual contributions of the 
diffusion/ convection, and pressure gradient terms in the momentum 
equation. In the step size study conducted here, 3 points in 
the layer were considered, namely the points: n = 6 which lies 

in the bottom of the viscous region, n = 8 near which the 
convection effects are negligible and the normal velocity V changes 
its sign, and n = 10 where all the terms effects are vanishingly 
small but in perfect balance. Step sizes of 0,4, 0.5, 0.6667 
and 1 . 0 are used to plot F at those three points . 

In Figures 12 and 13, the basically inviscid values of F(6) 
are presented. The three curves in Figure 12 are not straight, 
with the CC -scheme being- most nearly linear and, moreover, more 
accurate than the C-scheme. The UW-scheme gives better accuracy 
than the other two schemes , but is seen to be unreliable over 
this region of Ari . In Figure 13, both the RUDR and EB-schemes 
have very slight ciirvature with the RUDR scheme being the more 
accurate of the two. The TPE-scheme is found to give large 
errors and to lose its second order accuracy even with step size 


of 0.4. Figures 14 and 15 show the velocity at a viscous point, 
F(8) near the point where V changes its sign. Thus large 



41 


errors siioulS be expected to be -associated with the UW— scheme 
due to its required direction switching. This is the result 
observed in Figure 14 where again the CC“scheme is seen to - 
be the more accurate of the finite difference schemes. 

The exponential scheme results at the same viscous 
point n ~ 8 are shown, in Figure 15, in which it is seen that 
the EB-scheme has become more accurate than the others -and that 
the TPE-scheme encounters obvious curvature. Observations 
similar to those of Figure 12 for F{6) are obtained from Figure 
16 for P{10) . On the other hand. Figure 17 for the exponential 
scheme- results- for F'CIO)’- resembles . those -of Figure 15. for" F (-8)'. 

The F-velocity profiles for this case are presented in 
Tables 3 and 4 for step size of 1.0. For this moderate blowing 
case, oscillations are observed in the central difference 
calculations, not only near the outer edge as before, but also 
in the basically inviscid region close to the surface where 
again osciriatory: roo.tsy.to the difference equations exist. 

Errors' are observed in the calculations ' by the’-hotally exponential 
schemes- -(TPE and EB^schemes) in the region close ■'•to the” surface- 
where^^the-'Ve-locity profile is-linear.-.-i -These-' errors are- larger 
in the^case of' the: TPE: scheme*.. 

To this point, the following conclusions can be made based 
on the moderate blowing results presented here; 

1, The CC-scherae is more accurate than the C-scheme. 

2. Ap.art from the nonsystematic behavior of the UW^scheme 
due to the errors introduced at the region in which V changes its 
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sign, the UW-scheme gives reasonably accurate results. 

3. The RUDR-scheme represents a straightening of the 
CC- scheme. 

4. The RUDR and EB schemes are the two schemes that have the 
slightest deviation from the second order accuracy at la-rge step 
sizes. The EB scheme shows better accuracy as the injection 
velocity is increased, which should be expected since the transi- 
tion from the inviscid region created by the injection to the 
viscous region is better honored in the EB-scheme. 

5. The TPE-scheme is not comparable with ‘the other 
exponential schemes, neither in its accuracy nor ^in its behavior 
with the step size. The errors observed in the velocity profile 
in the inviscid region close to the surface are blamed for that 
and are assumed to affect the solution everywhere else. One 
other' possible reason for this limited accuracy could be that 
the integration of the continuity equation by the trapezoidal 
rule is not consistant v;ith the remainder of the method . 

6. The central difference solution is oscillatory whenever 
a condition ' similar to condition (16) is satisfied. Oscillations 
do not occur in any other scheme. However, obvious errors are 
observed when exponentials are used to approximate inviscid 
linear regions of the profiles . 

4 . Case III - Strong Blowing at an Axisymmetric Stagnation Point 

To examine the validity of the previous conclusions, a 
high injection rate = 25 is applied, with the pressure gradient 
■coefficient 6 = 0.5. As shown in .Figures 18 and 19, a large 
region (0 < n < 47) of inviscid flow is formed in which the 


velocity profile is essentially linear followed by a very small 
viscous region (47 < n < 53) and then a region Ct > 53) in 
which the' diffusion/ convection, and pressure effects are 
negligible. The points considered for detailed analysis of the 
numerical results are n = 48, 50, and 52 for the same reasons 
discussed in the previous case and the step sizes used are 
0.5, 0.8, 1.0, and 2.0. In Figure 20 the CC-scheme in the 
inviscid region shows a large turn and loss of its second order 
accurate nature for An > 0.8. The OW-scheme shows a nonsystematic 
behavior in that the slopes of its segments increase and decrease 
alternatively and the segment of smaller step sizes does not 
head towards the exact value of F(48). Large non second order 
truncation errors are shown for the C-scheme. 

In Figure 21, only the EB and the RUDR schems appear because 
the TPE scheme shovied very poor accuracy at that point in the 
profile. The EB and RUDR schemes give the same solution for the 
smaller step sizes and both deviate slightly from the second 
order result for An > 1* 

Figure 22 shows the finite difference results at the switch 
point (V - 0 at n - 50) , where it is seen that the UW-scheme 
gives again large truncation errors and nonsystematic behavior. 

The C and CC-scheraes are seen to lose their second order 
accuracy for An > 1. In Figure 23, the EE-scheme , as expected, 
shows better accuracy than the RUDR, while the TPE-scheme follows 
its usual behavior observed before. These observations are 
further demonstrated for the edge region (n = 52) in Figures 24 
and 25 where F(52) is presented. Large truncation errors and 
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overshoot (F > 1) are encountered with the C-scheme, while the 
CC-scheme shows a dramatic loss in second order accxiracy for 
large Ar). The UW-scheme is obviously nonsystematic ; the TPE- 
scheme does not behave in a second order fashion; the RUDR 
scheme is fairly accurate and it represents a straightening of 
the CC-scheme; and the EB-scheme is the most accurate and 
consistently maintains second order accuracy even for the large 
step size of 2.0, 

The velocity profiles presented in Tables 5 and 6 support 
the previous observations that the C-scheme is oscillatory in 
the inviscid region close to the surface and the region close to 
the outer edge whenever a condition similar to condition (16) 
is satisfied. Also, the errors observed in the linear region 
close to the surface in the TPE and EB calculations are observed 
again. These errors are due to approximating a straight line 
by exponential segments and can be reduced by refining the 
matching between these exponential segments and by considering 
the next term in the asymptotic expansions assumed for the 
coefficients a, b, and c in the differential equation (44) . 
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VII. CONCLUSION AHD RECOMMENDATIONS 

A study of numerical schemes for solving viscid/inviscid 
fluid flow problems has been presented. Truncation errors 
of the oscillatory type associated with the central difference 
solutions of these problems for large step sizes have been 
related to an explicit downstream contribution to the con- 
vection effects which violates the, physical process of con- 
vection. 

It has been fo\md that a natural correction to the central 
difference scheme (which is equivalent to an addition of 
artificial viscosity like term) has improved the accuracy of 
the central difference scheme, increased its range of possible 
applicability, and eliminated the n\imerical "wiggles”. A 
fourth order scheme obtained using this same line of thought 
has been recently presented by Wornom (Ref. 17). 

A second order accurate upwind scheme has been described 
and found to suffer a severe conservation problem at the 
switch point although it has shown reasonable accuracy away 
from that point. It is recommended that this point should 
be studied further to solve this switching problem. 

The exponential scheme proposed by Spalding and used by 
Roscoe has been interpreted in a way Csee Pruess., Hef. 8.) that 
has allov/ed a development of the exponential box -scheme. This 
interpretation gives also a rational, basis for introducing, higher 
order corrections. 

The numerical calculations carried out with the exponential 
box scheme have, been found to be very promising. Further study 
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is recommended both in the application of the scheme to a wider 
range of problems and in considering the higher order corrections 
that can be introduced to it. 
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APPENDIX A 


WIGGLES AND DIAGONAL DOMINANCE 

The numerical oscillations or "wiggles" associated with 
the central -difference calculations for large step sizes have 
been attributed, incorrectly, to round off errors produced due 
to the lack of diagonal dominance in the difference equation. 
In this appendix, it will be shown that the wiggles and the 
lack of diagonal dominance are two different problems that 
can occur independently. Consider the difference equation 
(8) which will be rewritten here for convenience 


a-| F. T + 
1 3+1 


a., P . + 
2 3 


“3 ^j-1 = “4 


(Al) 


where > 0 . 

The characteristic equation corresponding to equation 
(Al) is 

+ 02 A + = 0 (A2) 


whose roots are 


-02 + ^c^2 - 4a^a^ 
2a^ 


and 


~°2 ~ ~ ^“ l °3 

2ot 


(A3) 


(A4) 


Consider now the following cases : 



Case I 


“3 ^ ° 


In this case always negative while X^ is alv/ays 

positive. 

. 2 

Case II 02 >40^0^ > 0 and - C 2 > 0 ^ 

Both X^ and X 2 are negative- in this case .. 

2 

Case III 02 ^ “2 ^ ^ 

Both X^ and X 2 are positive,. 

2 

Case IV ©2 <4o^.o2 

Complex roots X^^ and X 2 are obtained. 

Thus oscillations are not expected in Case III, while they 
are expected in Cases I and II. In Case IV, X^ and X 2 can be 
put in the form pe“ and the exact solution for the difference 
equation will have the form 

Fj =• [R^ cosje + R 2 sinje] (A5)' 

Since the set of equations (Al) for j=l, ..., N+1 has real 
coefficients, R, and R., will be such that all the F.'s are 

'12 j 

real and the solution profile will have harmonic behavior. 
Diagonal Dominance 

Consider the system of equations (Al) for j=l, ..., N+1. 

The diagonal dominance condition for this system is 

I “2'! — l“il 


for all j . 
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This condition is the sufficient condition for the matrix 
of the coefficients to be nonsingular and hence invertible (Varga, 
Ref .13) .. Also, it has been shown by Keller (Ref. is.) ■ that 
condition. (A 6 ) is a sufficient condition for the Thomas 
algorithm to be nonsingular. 

It is clear from the previous analysis that diagonal 
dominance has to do with the magnitude of the coefficients 
° 1 ' “2 “3' wiggles have to do with their signs. 

E.g., if ct^ < 0 oscillations will take place according to 

Case I above, although diagonal dominance can still be satisfied. 

However, in the case when ci2 < 0 and > 0 , it can be 
shown that the diagonal dominance is a sufficient condition 
for the roots and real and, according to Case III, 

positive. 

Squaring (A 6 ) and noting that and are positive, 
one gets 

+ 20^00 

2 

(a^ - a^) ^ 0 

X^ and X2 are real. 


2 2 ^ 2 
0^2 > ai + a3 

so that 

“2 ' ^“l“3 - 
which proves that 
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APPENDIX B 


A GENERAL THREE POINT FINITE DIFFEBENCE EQUATION 


The approach used to obtain the midpoint upwind scheme 
in Chapter II can be used to produce a family of finite difference 
schemes if it is applied to a general point "s" a distance 
a Ari below point j where -1 ^ a +1. 

Consider equation (4b) and its derivative evaluated at 
point s. They give 


F - C F =0 
nn • n 
s 's 


(Bl) 


and 


F - C F =0 


(B2) 


For 0 £ a <_ 1, a Taylor's series expansion of F gives 


nn. 


F = F • . + aAri F +0 (An ) 

nHj ntg ntrig 

Substituting for F from (B2) one gets 

s 

(1 + a An C)F^^ = F + O(An^) 

ntg ntj 

Similarly, a Taylor's series expansion of F and the use 

^m 


(B3) 


of (Bl) gives 


[1 - a)AnC]F = F t O(An^) 

X n n 

s m 


(B4) 


Using central differences to represent F and F , the 
following second order accurate finite difference approximation 
of equation (4b) is obtained: 
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LJ- “ ariujr 


*r ann v- jr. 


+ [1 + i(l + 2a)AnC + aAn^C^3F._^ = 0 


(B5) 


A similar expression can be obtained for the case where 

-1 < a < 0 with F used instead of F . Combining both cases, 

— — n Ti 

n m 


one gets : 


[1 - |(l-2a)AnC - j(a - [ a | ) An ] Fj^^_ 

- [2 + 2aAnC + |a|An^C^]Fj 

+ [1 + j(l+2a)AnC + i(ct + I aj ) An^C^]Pj_^ = 0 (B6) 


for -1 _< a £ 1. 

The corresponding second order error is given by 

" = -'ll ^nnn ^'l4-^l'|- (B7) 


where F and F are evaluated at some point and can be 

nnn nnnn 

related to each other through the differential equation (4b). 
so that 

e = [-(^ + I .2) + ly + 1 (| - |a|)2]An" (B8) 

Note that when a = 0, equation (B5) reduces to the central 
difference equation (7) , when ct “ ^ the branch of the midpoint 
upwind scheme corresponding to C > 0 is obtained and, when 
a = - — the branch corresponding to C < 0 results in. 
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APPENDIX C 


ERROR ANALYSIS FOR EXPONENTIAL SCHEMES 


1. Error Analysis of Equation (46) 

-a . Ati/ 2 -b.An/2 

Multiplying equation C46) by e ^ e ^ results 

in the following equation 

-a.An/2 -b.Ari/2 a.Ari/2 “b.An/2 -a. An/2 b.An/2 

e ^ e ^ 5*..^ - [e ^ +e^ e^ ]F. 


a.An/2 b.An/2 c. -a.An/2 -b.An/2 

+ e 3 g 3 p = e 

JO 

a.An/2 -b.An/2 -a.An/2 b.An/2 a.An/2 b.An/2 

-[e^ +e^ ]+e^ 


Expanding the exponentials for small An one gets 


-a.An/2 -b.An/2 , 

i ^ e ^ = 1- -x- An(a..+b.)+ ~ An^ (a .+b . ) An"^ (a .+b . ) “”+0 (An *) 


12 2 13 3 

‘ "^'i.+b.) - ~ An (a.+b.' 

J J 48 '^jj- 


3.. An/2 b.An/2 , 12 213 34 

= 1+ i An(a.+b.)+ ^ An (a.+b.)^+ ^ An (a.+b.)^+0(An ) 
2 3 3° 33^° 33 


-a .An/2 b .An/2 
e ^ e ^ 


= 1- y An(a.-b,)+ ^ An^(a.-b.)^- ^ An^ (a .-b . ) ^+0 (An^) 

2 Djy J J ■ 48 J J 


a.An/2 -b.An/2 


= 1+ J An(aj-bj)+ ^ An^(a^-bj)^+ An ^ (a^-b^ ) ^+0 ( An^) 


■ (C2) 


Substitution in the right hand side of equation (Cl) gives 


R.H.S. = An^ Cj + 0 (An'^) 


(C3) 


Expanding and ^j_j_ about point j one gets 

= F. + AnF + i An^ F + An^ F + OCAn^i 

3+1 J ' nj 2 ' ntj 6 nnnj 


rssss 


S3 



F. , = F . - AnF^ + j A.ti^ F - i An^ F + O(An^) 

3 -*- J rij ^ nrij o rmrij 

Siibstitution in the left hand side of equation (Cl) gives 
coeff. of F^ = 1 - I An.(aj+bj)+ | ati^ (a^+bj) An^(aj-f-bj 

- 2 - J An^(a.-bj)^ 

+ 1 + ^ Ar>(aj+bj)+ ^ At]^ (a^+bj ) ^ An^(aj+bj 

= An^ a^bj 

coeff. of F = Ar\“ y ATi^(a.+b.)+ i An^(a.+b.)^ +0 (An ') 

^ J J o J 3 

- An- j An^(aj+bj)- An^(ajtbj)^ +0(An'^) 

2 4 

= “ An (a^-i-bj) +0(An') 

coeff. of F^^_ =1 An^ - I An^ (a^+b^) + 0(An^') 

+ y An^ + An^ (a^+bj) + 0(An"^) 

2 4 

= An + 0 (An *) 

coeff. of F = y An^ + O(An^) 

nnn j ° 

- I An^ + O(An^) 

= 0 (An *] 


(C4) 

) ^+0 ( An ’) 

+0 (An^) 
) ^+0 (An') 
+0 ( An^) 
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so •cna'c 


L.H.S = An^[F -Ca.+b.)F + a .b .F .]+ 0 (An^) 

TiTlj 3 D Uj 3 3 3 


(C5) 


Equating (C3) and (C5 ) , one gets 

F - {a.+b.)F + a.b.F. = c. + 0 (Ar)^) 

ntj 3 3 Hj 3 3 3 3 

which shows that equation (46) is a second order accurate 
representation of equation (44) . 

2 . Error Analysis of Equation (50) 

-b . An/2 

Multiplying equation (50) by e -* results in 


-b . An/2 b . An/2 

-e ^ F . + e ^ 

3 


j-1 a .b . 

3 3 


c . b . A n / 2 — b . A n / 2 

(e 3 - e ^ ) 


(C7) 


Expanding the exponentials in the right hand side, one gets 

c 

R.H.S. = — {[1+ ~ Anb. + i An^ b? + O(An^)] 

ajOj 3 B j 

- [1 - i Anbj + I An^ b>j + 0 (An^) ] } ^ 

c . ^ 

= An — ^ — H 0 ( An ) 

3, > 

3 


(C8) 


.Expanding the exponentials in the left hand side and P-_-, about 

3 ~ J. 

point j one gets 


obigwalpagb® 
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L.H.S. = -[1- i An b . + i An^ + 0(An^)]F. 

^ J O J J 

+ [1+ ^ Anb.+- An^bJ+0(An^) J [F.-AnF + ^ An^F +0(An^)3 
^ 3 ° 3 3 Tij ^ nHj 

= Anb.F.-An(l+ i Anb.)F + ^ An^F + O(An^) (C9) 

3 3 3 Hj ^ nnj 

So that, upon equating CC8) and (C9) there follows; 

c - , « 

-P^ + b.F. = ^ i An(F^„ - b . F ) + 0{An ) (CIO) 

rij 3 3 a^ 2 nnj 3 ^ij 

which shows that equation (50) is a first order accurate repre- 
sentation of equation (49) . It is second order accurate only if 
the coefficients a, b, and c are constants. 

3. Error Analysis of Equation (58) 

.An error analysis for equation (58) has been carried out by 

expanding its coefficients and reducing it to the three point 

difference equation obtained by applying the- Keiler'-s box scheme 

to the same differential 'equation (44) , The difference 

2 

between the two schemes rs found to be 0 (An ) which shows that 
the EB scheme is second order accurate since the Keller's box 
scheme is known to be second order accurate. 

4 . Error Analysis of Equation (59a) 

-bj^An/2 

Multiplying by e , one gets 

-b.,An/2 ,b, An/2 c, b, An/2 -b., An/2 

-e . F. + e F. , = — ~ [e - e ] 

3 3-1 aj_b^ 
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(Cll) 



R.H.S. 




{ [1+ J An + I An^ + 0 (An.^) ] 
[jL - ^ An + |- An^ b^ + OCAn^)] 


L.H.S. 


An + 0 (An^) 

ai 

-[1- I Anb^+ I An^bJ +0 (An^) ] EF^+ 

+-[1+ I Anb^+ I An^b^ +0(An^}] [F^^- 
AnbiF^ - AnF^ + O(An^) 


(C12) 

- AnF„ + i An^F^^ +0 (An^) ] 

i AnF^ + J- An^F +0 (An^) ] 

z x\ o Ti 

'm ' m 

(C13) 


Equating (Cl2) and (Cl 3) , one gets 

2 

-P + b,F = “ + 0(An ) (C14) 

T| JL Xtt ^ ^ 

m 1 

which shows that equation (59a) is a second order accurate 
representation of equation (49) , 
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APPENDIX D 


A SECIOND GENERATION EXPONENTIAL BOX SCHEME 


Let the coefficients a, b, and c in equation (44) be 
approximated to first order by their values at the grid point 
j and the solution (45) be valid in the interval between the 
two midpoints m and n. For convenience, equation (45) will 
be rewritten in the form 
a^n 

F = A 2 e + B 2 e- + ^2 ' 


where 


^2 a«b 




(D2) 


and the sxobscript "2" refers to point j. 

Similar expressions will be valid at midpoint intervals 
below and above the previous interval with the coefficients 
a, b, and c approximated by their values at points j.-L and 
j+1, respectively. 


Equation (Dl) can be used to evaluate F., P , F , F 

J m n' n 

and F in the forms 


m 




^2 

+' 

®2 

+ 

^2 

(D3) 

''I 

’ H 



-1 







+ 

B2 t 2 ^ 

+ 

^2 

(D4) 



-1 






^n 


^2 ^2 

+ 

^2 

+ 

^2 

(D5) 

F 

= 

^ 2^2 ^2 


^ 2 ^ 2"^2 



(D6) 


^Tv ^2'^2^2 ^2®2^2 
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where r 2 


-a-^Ari/2 b„An/2 

e and t 2 = e 


(D8) 


Three expressions for F . , , F , and F can also' be obtained 
using the expression for F in the lower interval. These 
expressions will be similar to the R.H.S. of'(D3) , (D5) , and (D7).. 

with the subscript "2" replaced by "1" to refer to point j-1. 

I 

On. the other hand the expression for F in the upper interval 

will give for and F^^ expressions that are similar to 

the R.H.S. of (d3 ) , (d4 ) , and (D6) if the subscript "2" is replaced 

by "3" to refer to the point j+1. 

The resulting set of 11 equations is sufficient to eliminate 

A, , B, , A„, B-,, A_, F , F , F , and F . and an .equation 

1 ± z z z j m n 

that involves ^ j ^j+1^ only is obtained 

This resulting equation is tiie difference equation for 
this scheme. It has the form 

^^l“^l^®3^jtl j + ^^3“^3^®1^ j-1 

= (g^-h^)d2 - (g^h^-g^h^) f2 + (h2-g3)d^ (D9) 


where 

di = (a^-b3_r3_t^)t2f2 + [{a^-b^)t^ - (a^-b3_r^t3_) ] t2f;j_ 

d3 = (a3r3t3-b3)r2f2 + [(a3-b3>r3 - (d3r3t3-b3) ] r2f 3 

e^ = t^t2(a3^-b^) 

®3 ^ 2 ^ 3 ^^ 3 "*^ 3 ^ 


ORIGINAL PAGE IS 
OP POOR QUALITY 


61 



9l = [(a^-b^r^t^) - a2 (1 - 

53 = [(a3r3t3-b3> +33(1 - r3tj) ] 

**1 “ ■ ’=■ 2 '^ - 

= [(a^r^t^-b^) + b2 (1 - r2t2)]r2t2 (DIO) 


The scheme given by equation (D9) has the advantage over 
the scheme given by equation (58) that no values need be eval- 
uated at the midpoints m and n. 

Instead of using the Thomas algorithm to invert the .set 
of three point difference equations given in (D9) a simpler 
inversion scheme can be designed to first evaluate the coef- 
ficients Aj and (j=l, N+1) and then the values of 

F and its derivatives can be calculated directly from (Dl) in 
the appropriate interval. Such an- algorithm eliminates the 
tedious algebraic steps necessary to derive the tridiagonal 

equation (D9) and is described here for future reference. 

/ 

The requirement that and F^^ be continuous at point ra, 

the midpoint between the j-1 and jth points, provides two 


equations that involve A . , , B . , , A . , and B . that can be solved 

^ I -1 ' ■] -1 ' n ' T 


for A . and B . . Thus one can write 
3 5 . ■ 


A. 

3 


“ij ^j-i ^ ®i-i 



“ 2 j ^j -1 + 62 'j ®j-l + ^ 2 j 


(Dll) 

(D12) 


where the a’s, 6 's, and y’s are all known. 
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Now 
the form 


assvuning that Aj and (j=l 


,N+1) can be written in 


A • ” |i T • At 

3 "^13 1 

"Ij 

(D13) 

Bj = pj, + 

"2j 

(D14) 


and substituting ^j-1 (Dll) and (D12) one gets 


= (a..u 




B. 

3 


= (a 




)A, 


<“2j'’l,j-l'"®2j'’2,j-l+’f2j> 


Comparing (D15) and (D16) with (D13) and (D14) the following 
recurrance relations for a.nd obtained 


^Ij ®lj ^l,j-l ^Ij ^2,j-l 

^Ij = “Ij ^l,j-l »lj ^2,j-l + ^Ij 

^2j " “2j ^l,j-l ^2j ’^2,j-l 

v„ . =s . Vt . , + . v_ . , + . 

2j 2j l/:-l 2j 2,3-1 2j 

At j = 1 the following relations are known 

“i.i “ ^ 

^>1 1 = 0 

' X • 

ORIGINAL PAGE B’ ' 
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(D17) 


(D18) 

(D19) 

(D20) 



while at j = N+1 the following relations are valid 


^+1 

^1,N+1 

A^ + 

^1,N+1 

(D21) 

®N+1 

^2, N+1 

Ai + 

'*2, N+1 

(D22) 

^N+1 

■^+1 

®N+1 

^N+1 

(D23) 


In (D20) and (D23) , and are assumed to be known as 

boundary conditions of the differential equation under consideration. 
Other types of bounda3ry conditions can be treated using expression 
(Dl) for F to give two relations involving (A^, B^) and 
to replace (D20) and (D23) respectively. 

The calculation procedure will be as follows : 

1. Use (D20) to obtain ^2 and by comparing it with 

(Dl4) for j=l- 


2. Use the recurrence relations (D17) to calculate y. 


Ij' 


^ij , U 2 j/ N+1 knowing their values for 

j=l. 


3. Knowing the values of ''l,N+l' “2,N+1' 

solve equations (D21) , (D22) and (D23.) to determine . 

4. Calculate the ^j's B^'s, j~2 , N+1 using (D13) 

and (D14) . 

5. Calculate F and its derivatives using (Dl) . 
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n 

C 

CC 

UW 

0,0000 

0,0000000000000 

0. ooooooooooooo 

0 . ooooooooooooo 

0 , 5000 

0,2357019175320 

0.2330208993764 

0,235107681 io25 

1 .0000 

0, 4645601396245 

0.45935227761,07 

0,4646173272125 

1 ,5000 

0.6681224 JB0019 

0.6609601688855 

0,6682772996516 

2,0000 

0,8250732563570 

0,8170778673049 

0.8232959966979 

2.5000 

0,9248680267167 

0,9175782221355 

0.9205713905362 

3 ,0000 

0.9749164338430 

0.9696712237985 

0,9700708297525 

3,5000 

0.9938657448810 

0,9910125484759 

0,9905255765472 

4,0000 

0. 9989977732033 

0,9978690991754 

0,9974650685493 

4,5000 

0. 9999093638850 

0,9995958195990 

0.9994204753158 

5,0000 

0.9999977293927 

0.9999382939866 

0,9998ftb06o9o93 

5.5000 

1.0000000765307 

0.9999923165958 

0.Q99980555/|q<);a 

6,0000 

0 ,9999q99934o52 

0,9999992026792 

0.9999<?7o966S43 

6,5000 

1,0000000008841 

0,9999999289259 

,0.9909996179303 

7,0000 

0.9999699998435 

0,9999999943462 

0.0999900553871 

7 ,5000 

1 . 0 OOO 0 OOOOO 34 O 

0,9999999995812 

0.9999999953493 

8,0000 

0,9999999999915 

0,9999999999699 

0.9909999995647 

8,5000 

1 ,0000000000025 

0.9999999999978 

0,9999999999632 

9,0000 

0,9999999999993 

0.9999999999998 

0,9909999999972 

9,5000 

1 ,O0DOo0O0o0OO4 

1, ooooooooooooo 

0.9999909999998 

t 0,0000 

1 .ooooooooooooo 

1 , OOOOOOOOOOOOO 

1. OOOOOOOOOOOOO 


-TABLE 1: BLASIUS FLOW 

VELOCITY PROFILES FOR FINITE DIFFERENCE SCHEMES 
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EUDR & TPE 

EB 

0.0 

0.0 

0,0 

0.5000 

0,2330545736248 

0.2354540274489 

1,0000 

0,4594177063782 

0,4624499528736 

1.5000 

0.6610505304537 

0,6632486160505 

2,0000 

0.8171816649639 

0,Bl7'»955i91227 

2.5000 

0,9176822625149 

0.9175423945048 

3.0000 

0.9697617889022 

0,9693773037453 

3.5000 

0,9910775955934 

0,9908197225985 

4.0000 

0,9979055192255 

0.9977999677228 

4,5000 

0,9996112085620 

0.9995811491306 

5.0000 

0,9999431583436 

0,9999369107351 

5,5000 

0,9999934731040 

0.99999250410‘52 

6,0000 

0.9999994125816 

0,9999992990353 

6.5000 

0,9999999586259 

0,9999999484074 

7,0000 

0,9999999977220 

0,9999999970308 

7,5000 

0,9999999999020 

0.9999999998658 

8.0000 

0.9999999999967 

0,9999999999955 

8,5000 

0 ,9999999999999 

0,9999999909999 

9.0000 

1 .OOOOOOOOOOOOO 

1 , oOOOfiOOOoOnOO 

9,5000 

1 .OOOOOOOOOOOOO 

1 .OOOOOOOOOOOOO 

10,0000 

1 .OOOOOOOOOOOOO 

1 .OOOOOOOOOOOOO 


TABLE 2: BLASIUS FLOW 

VELOCITY PROFILES FOR THE EXPONENTIAL SCHEMES 
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n.noooooooouooo 
0, 1 25000 I 385 I 62 
0,?/J999970509i46 
0.3750012121533 
0.a99'?9i|957«509 
0, 6250309063076 
0.7^95805282511 
0.81555708850288 
0.9511730093339 
0.9912926329808 
0.9999098573588 
1 ,09OO) 16855790 
0,999997l7278ii0 
1.0000009069873 
0.9999996238077 
1 , 000000l8i!8fc76 
0.999999916‘1921 
1 .0000000599589 
0,9999o997997«7 
1.0000000297870 

i.ooooooooooooo 


0, ooooooooooooo 
0, 1 209999975309 
0,2il999996£)5o76 
0,37«999«929750 
0,il99992348929« 
0,6248863622985 
0,7466j80257o22 
0 , 86H009C 195098 
0,95olo65bo5221 
0.9898745907980 
0.9959944930587 
0 , 9999004691035 
0 , 9999963.150787 
0,9999996758327 
0,9999999594763 
0,9999999932714 
0.9999999986084 
0,999999999661 7 
0.9999P99999103 
0.9999999999765 

1 .OOOOOOOOOOOOO 


0 , OOOOOOOOOOOOO 
0, 1249999909974 
0.2499998742456 
0.3749984090375 
0.4999814704162 
0.6248108847319 
0,7484101762561 
0.8647799101187 
0,9539534103005 
0,9907807947510 
0,9988564309871 
0,9999o59459o 89 
0.9999945558*09 
0,9999997674294 
0,p999p999?3875 
0.9999999998033 
0,9999999999959 
0.9999999999099 

1.0000000000000 

1 .0000000000000 

1 .OOOOOOOOOOOOO 


TABLE 3 

MODERATE INJECTION AT AN AXIS^TMMETRIC STAGNATION POINT 
VELOCITY PROFILES FOR FINITE DIFFERENCE SCHEMES 


0.0 

0,1249999999508 
0,2499999973900 
0,3749998831525 
0,4999960124792 
0,6249072666931 
0,7486703129367 
0,8640974059548 
0.9501418039624 
0.0699214551605 
0,9990551540856 
0,9999627371932 
0,99999t)al8b416 
0.0999999965059 
0,9999999999920 

1 .OOOOOOOOOOOOO 

1 .0000000000000 

1.0000000000000 
1 .ooooooooooooo 

1 .0000000000000 

1 .ooooooooooooo 


0,0 

0.1255323U9837 

0,2515652915513 

0.3780471033645 

0.5048612680820 

0.6316647394431 

0.7567270901976 

0.8715107771055 

0 9544158150455 
0.9911497573511 
0,9992034627362 
0.9999699736813 

0.9999995509957 

ft. 9999099974 1 14 

0) 9999999999943 

1 ) OOOOOOOOOOOOO 

1.0000000000000 

1 .OOOOOOOOOOOOO 

1.0000000000000 

1 .ooooooooooooo 

1.0000000000000 


0.0 

0,1251331647276 

0 2502705990423 

0)3754i52395l57 

0 5005658740969 
0.6256233708841 
0.7493524451070 
0,8641102486477 
0,9492358621257 
0.9892008319580 
o)o988872841214 
0.9‘’99502945726 
0.0999991012758 

0.9999999937433 

0)9999999999838 

I.OOOOOOOOOOOOO 

I.OOOOOOOOOOOOO 

1 .OOOOOOOOOOOOO 
I , oOOOOOOOoOOOO 
1 .ooooooooooooo 
I .ooooooooooooo 


TABLE 4 

MODERATE INJECTION AT AN AXISYMMETRIC STAGNATION POINT 
VELOCITY PROFILES FOR THE EXPONENTIAL SCHEMES 
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c 


cc 


uw 


Q. 01)00 Q.OOOOOOOOOOOOO 0 , 0 0 0 0 0 0 0 00 0 00 0 O.oOOOoOOOOOOOO 

2.0000 O.o399oitbi91 128 0 . o3999999852&6 0 ,n«00000000000 

i^.OOOO 0,0800000871 197 0,0799999968522 0 .0800000000000 

6.0000 0,1 199900901971 0 . \ ) 99999942626 0 , t SOOOOOOOOOOO 

8.0000 0,1800010729529 0.5 599999912057 0 . 5 800QOOO O'OOOO 

10.0000 0.1999933981 039 0. 5 ^99999872766 0 ,2000o0000‘0000 

12.0000 0,2000018778280 0,239999982196 5 0 , pOOOnoOOOOOOO 

14.0000 0,2799924667368 0,2799999755761 0,2800000000000 

16.0000 0,3200o2'^7921 1 1 0.3199999668305 0.3200000000000 

18.0000 0.3599911701 180 0.3599999553297 0 .36000000 OOOOO 

20.0000 0,4000045428948 0,3999999397991 0.4000000000000 

22.0000 0.4399892873270 0,4399999186 5 33 0.44001)00000000 

24.0000 0,4800068696053 0,4799998892170 0 . 4800000000000 

26.0000 0,5199864018206 0.5199998476014 0.5200000000000 

28.0000 0.5600105555772 0,5599997872657 0,5600000000000 

30.0000 0.5999816505634 0 , 5‘7999969 725 1 3 0 .6000000000000 

32.0000 0.6900168896633 0.6399995582156 0,6400000000000 

34.0000 0.6799730012435 0.6799993340641 0,6800000000000 

36.0000 0,7200292115698 0,7199989526665 0.72OOO00000000 

38.0000 o!75995U697196to 0,7599983568631 0 .7600000000000 

40.0000 o|80005822?5960 0 , 7'799968o34 1 83 0,7*799999999994 

42.0000 0.8399|3«88l2l50 0.83999 368339 30 0,8399999998945 

44.0000 0|8801548127682 0 .8799gii860 1666 0.87999^9826392 

46.0000 0,9196792203024 0 , 9 1 99s 1 07 1 9 q 1 4 0.9 5 99980557216 

48.0000 0,9608975687057 0,9597100879697 0.9593707621555 

50, .0000 0^9935820005289 0.99 31702198606 0 ,995892 1 8o993 1 

52.0000 1.0008980904585 0.9997149261997 0.9998778987386 

54.0000 0.9996811278535 0,9999503497977 0,9999^83574438 

56.0000 1,0001560202453 0 , 99998400o509fa 0 ,'’99909987 19 1.2 

58.0000 0,9*799066937334 0.9999930172200 0.9999999999561 

60.0000 1.0000595721305 0,9999963653360 0,9999999999998 

62.0000 0.9999563522191 0.9999978724550 I , qOOOOOOOOOOOO 

64.0000 1,00003055o 5276 0.9999986521368 1 . oOOOoOOOOOOOO 

66.0000 0,9999744813318 0.9999990946317 1 , OOO OOOOOOOOOO 

68.0000 1 . OOOO 1 8 I5l 6 1 1 0 0,9999993644662 1 , oOOOoOOOOOOOO 

70.0000 0,9999829476229 0.9999995385277 1 .qOOOOOOOOOOOO 

72.0000 1.00001 17006454 0,9999996560146 5 .qO OOOOOOOOOOO 

74.0000 0.9999875111227 0,999999738-3396 1 , QOOOOOOOOOOOO 

76.0000 1.0000o78837714 0.9999997978756 1 .OOOOOOOOOOOOO 

78.0000 0,9999902057064 0.9999998421129 I .OOOOOOOOOOOOO 

80.0000 1.0000054074186 0,9999993757682 1 .qOOOOOOOOOOOO 

82.0000 0,9999918952477 0.9999999059527 1. OOOOOOOOOOOOO 

84.0000 1,0000036835870 0.9999999226055 I . nOOOoOOOOOOOO 

86.0000 0,9999929967801 0,9999999392623 1 , oOOOOOOOOOOOO 

88.0000 1,0000024141944 0*9999999528787 t , qOO OOOOOQOOOO 

90.0000 0,9999937515261 q ‘ g999999o4 1689 i . qOOOOOOOOOOOO 

92.0000 1,0000014351876 0*9999999736543 1 , oOOOoOOOOOOOO 

94.0000 0,9999942256535 0,9999999317216 1 , qOOOOOOOOOOOO 

96.0000 1 .0000006502884 0 1 909999988662 1 1 . oOOOqoOOOOOOO 

98.0000 0,9999945553756 o ,9999999940978 1 . qOOOOOOOOOOOO 

100.0000 1 .OOOOoOOOOOOOO 1 , o 0 0 0 Q 0 0 0 0 0 0 0 0 I , 0 0 0 0 QQO 0 0000 0 

TABLE 5 

STRONG BLOWING AT AN AXISl&IMETRIC STAGNATION POINT 
VELOCITY PROFILES FOR FINITE DIFFERENCE SCHEMES 
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n 

RUDR 

TPB 

ER 

0.0 

2,0000 
«,0000 
Cl. 0000 

fl.oooo 
10,0000 
12,0000 
la.oooo 
16,0000 
18,0000 
20 .0000 
22.0000 
20.0000 
26.0000 
28.0000 
30.000C 

32.0000 

30.0000 

36.0000 

38.0000 
00.0000 
a2,0000 
^0,0000 

46.0000 

48.0000 

50.0000 
• 52,0000 

54.0000 

56.0000 

58.0000 

60.0000 
62.0000 

64.0000 

66.0000 
68,0000 

70.0000 

72.0000 

74.0000 

76.0000 

78.0000 

80.0000 
82.0000 

84.0000 

86.0000 

58.0000 

90.0000 

92.0000 

94.0000 

96.0000 

98.0000 
100,0000 

0,0 

0,0400000000000 
O.O^OOOOOOOOOOO 
0.1200000000000 
0, 1 OOOOOOOOOOOO 
0.2000000000000 
0 . 2^00000000000 
0.2800000000000 

0.3200000000000 

0.3800000000000 

0.4000000000000 

0.4400000000000 

0.4600000000000 

0,5200000000000 

0.5600000000000 

0,6000000000000 

0.6400000000000 

0,6800000000000 

0,7200000000000 

0.7600000000000 

0.8000000000000 

0,8400000000000 

0 .ft7q99999<59g99 
0.9199999785978 
0.9599317532875 
0,99325221 1 3264 
0,9999367431587 
0.9999999854915 
0 ,9999099909999 
1 . OOOOOOOQOOOOO 

1.0000000000000 

1.0000000000000 

1 .0000000000000 

1 .0000000000000 

1.0000000000000 
1 ,0000000000000 

1.0000000000000 

1 .0000000000000 

1 .0000000000000 

1 .0000000000000 

1 .OOOOOOOOOOOOO 
1 .OOOOOOOOOOOOO 
1 .OOOOOOOOOOOOO 
1 .OOOOOOOOOOOOO 
1 .OOOOOOOOOOOOO 
1 .OOOOOOOOOOOOO 
1 .OOOOOOOOOOOOO 
1 .OOOOOOOOOOOOO 
I .OOOOOOOOOOOOO 
1 .OOOOOOOOOOOOO 
1 .OOOOOOOOOOOOO 

0.0 

0.040030756o792 

0.o8009256o8737 

0 , 1 20 1 849905026 

0.160io89l'27o62 

0 .20046448726o5 

0.2406521713189 

0.2808725254664 

0,321 1262212638 . 

0,3614140500788 

0.4017369331309 

0.4420959321954 

0.4824922597755 

0.5224272861505 

0,5634025378059 

0.6039196755681 

0,6444804270795 

0.6850864163865 

0.7257387541931 

0, 7664370380782 

0.8071767653865 

0,8479419049149 

0,88867978(16307 

0.«291 924545652 

0.9683961219803 

0.9960251114147 

0 , 9999742948492 

0 . 9999999959049 

J .OOOOOOOOOOOOO 

1.0000000000000 

1 .0000000000000 

1.0000000000000 
l.ooooooooooooo 

1.0000000000000 

1 .OOOOOOOOOOOOO 

1 .0000000000000 
1 .OOOOOOOOOOOOO 
1 .OOOOOOOOOOOOO 

1 .0000000000000 
t .OOOOOOOOOOOOO 
1 .OOOOOOOOOOOOO 

1 .0000000000000 

1 .0000000000000 
1 .OOOOOOOOOOOOO 
t .OOOOoOOuOOOOO 
1 .OOOOOOOOOOOOO 
1 .OOOOOOOOOOOOO 
1 .OOOOOOOOOOOOO 
1 .OOOOOOOOOOOOO 
1 .OOOOOOOOOOOOO 

1.0000000000000 

0.0 

0 . 0^00006167749 

0.0800012187300 

0,1200017875510 

0.1600023039530 

0.2000027470678 

0.2400030937425 

0,2800033177037 

0.3200033885324 

0.3600032703770 

0 .4000029203098 

0.4400022862057 

0.4800013039912 

0,5199998940945 

0.5599o;9S696o5 

0,S999q53o77714 

0.639991971 1466 

0.6799675794772 

0.7l998l98609o5 

0.7599750306494 

0.7999668431905 

0,8399^37283480 

0,87995o48l4003 

0.9199320069041 

0.9599439967085 

0.99232197 15474 

0,999B746«22664 

0,9999999902784 

0 . 6999999999999 

1 .0000000000000 
1 .OOOOOOOOOOOOO 
1 .OOOOOOOOOOOOO 
1 .OOOOOOOOOOOOO 

1.0000000000000 

1 . OOOOOOOOOOOOO 
t .oOOOoOOOOOOOO 

1.0000000000000 

1 .OOOOOOOOOOOOO 
t .OOOOOOOOOOOOO 
1 .OOOOOOOOOOOOO 

1 .0000000000000 
1 .OOOOOOOOOOOOO 
1 .OOOOOOOOOOOOO 

1 .0000000000000 

1 .0000000000000 

1 .0000000000000 
1 .OOOOOOOOOOOOO 
1 .OOOOOOOOOOOOO 
1 .OOOOOOOOOOOOO 
1 .OOOOOOOOOOOOO 
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FIGURE 1. VISCID/INVISCID FLUID FLOW 
OVER K REENTRY BODY 
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FIGURE 3: CHARACTERISTIC ROOT OF FINITE DIFFERENCE • EQUATIONS 







FIGURE 5 : B LAS I US PLOW-MOMENTUM BALANCE 
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FIGURE 7 


BLASIUS FLOW-STEP SIZE STUDY FOR THE EXPONENTIAL 
SCHEMES, F AT n = 1.0 
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FIGURE 9 


BLASIUS PLOW-STEP SIZE STUDY FOR THE EXPONENTIAL 
SCHEMES, F AT n = 2.0 
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FIGURE 11: MODERATE INJECTION AT AN AXISYMMETRIC STAGNATION 

POINT-MOMENTUM BALANCE 
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FIGURE 12: " "MODERATE ' INJECTION- 'AT AN AXISYMMETRIC STAGNATION 
POINT-STEP SIZE STUDY FOR FINITE DIFFERENCE 
SCHEMES^ F AT n = 6.0 










FIGURE 14; MODERATE INJECTIOM AT AN AXISYMMETRIG STAGNATION POINT-STEP 
. SIZE STUDY FOR FINITE DIFFERENCE SCHEMES, F AT n = 8.0 
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FIGURE 15: MODERATE INJECTION AT AN AXISYMMETRIC STAGNATION POINT-STEP SIZE 

STUDY FOR THE EXPONENTIAL SCHEMES, F AT n = 8.0 
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FIGURE 17; MODERATE INJECTION AT AN AXISYMMETRIC STAGNATION ROINT- 
STSP SIZE STUDY FOR THE EXPONENTIAL SCHEMES, P AT n = 10. 
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FIGURE 18: STRONG BLOWING AT MI AX 

POINT-VELOCITY PROFILE 
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FIGURE 20; STRONG ’BLOWING AT AN AXISYMMETRIC STAGNATION POINT-STEP SI2E 
STUDY FOR FINITE DIFFERENCE SCHEMES, F AT ri = 46 






FIGURE 21: STRONG BLOWING AT AN AXISYMMETRIC STAGNATION POINT-SOTEP SIZE 

STUDY FOR THE EXPONENTIAL SCHEMES, P AT n = 48 










FIGURE 24; STRONG BLOWING AT 2\n AXISYMMETRIC STAGNATION POINT-STEP SIZE STUDY 
FOR FINITE DIFFERENCE SCHEMES, F AT n = 52 






FIGURE 25: STRONG BLOWING AT AN AXISYMMETRIC STAGNATION POINT-STEP SIZE STUDY 

FOR THE EXPONENTIAL SCHEMES, F AT n * 52 




